import os
import sys
import numpy
import math
import osr
import gdal
from gdalconst import *
import psycopg2
import psycopg2.extras
import datetime
from datetime import timedelta
def quadratic(aod,temp,avg_temp,avg_rh,avg_preci):
	pm25 = 14.5824131094822*pow(aod,2) + 0.00717672042227793*pow(temp,2) + 0.0436745050768178*pow(avg_temp,2) + 0.110567956642307*pow(avg_rh,2) + 0.0223066644364859*pow(avg_preci,2) + 0.204032026317446*aod*temp - 3.64981813542097*aod*avg_temp - 1.72577215701688*aod*avg_rh + 2.00804343652011*aod*avg_preci + 0.0700766661716725*temp*avg_temp + 0.0488985684555605*temp*avg_rh + 0.0880021896381959*temp*avg_preci + 0.0243399434140353*avg_temp*avg_rh - 0.0779563259677111*avg_temp*avg_preci - 0.192243725097092*avg_rh*avg_preci + 166.100857354171*aod - 10.1354884677253*temp - 25.1010391955021*avg_temp - 29.3493345538898*avg_rh - 11.3435053380987*avg_preci + 2923.076893329164
	return pm25

def get_year(filename):
	previous = 0
	count = 0
	year = "-1"
	for i in range(0, len(filename)):
		if filename[i] == '.':
			count += 1
			if count == 2:
				year = filename[previous+2:previous+6]
				break

			previous = i;

	return year;

def main():
	con = None
	try:
		con = psycopg2.connect(host='192.168.3.4', database='fimo', user='rasdaman', password='rasdaman') 
		cur = con.cursor(cursor_factory=psycopg2.extras.DictCursor)
		
		cur.execute("SELECT mod04.aqstime, mod04.filename,mod04.filepath, mod07temp.filename as temp_filename,mod07temp.filepath as temp_path from (SELECT aqstime, filename,filepath FROM apom.satresampmyd04 where aqstime < '2014-09-12 00:00:00'::timestamp) as mod04 inner join (SELECT aqstime, filename,filepath FROM apom.satresampmyd07temperature) as mod07temp ON (mod04.aqstime = mod07temp.aqstime)")
		rows = cur.fetchall()
		for row in rows:
			aqstime = datetime.datetime.strptime(str(row["aqstime"]), "%Y-%m-%d %H:%M:%S") + timedelta(hours=7)
			month= aqstime.month
			year = aqstime.year
			aod_filename = row["filename"].strip()
			aod_path=row["filepath"].strip()
			aod_path=os.path.join("/var/www/html",aod_path)
			aod_mask_file=aod_filename[:-4]+"_mask.tif"
			aod_mask_file=os.path.join(aod_path, aod_mask_file)
			
			temp_filename = row["temp_filename"].strip()
			temp_path = row["temp_path"].strip()
			temp_path = os.path.join("/var/www/html",temp_path)
			temp_mask_file=temp_filename + "_T_10km_mask.tif"
			temp_mask_file=os.path.join(temp_path, temp_mask_file)
			
			aod_dataset = gdal.Open(aod_mask_file)
			aod_data = aod_dataset.GetRasterBand(1).ReadAsArray()
			
			temp_dataset = gdal.Open(temp_mask_file, GA_ReadOnly)
			temp_data = temp_dataset.GetRasterBand(1).ReadAsArray()

			cols = aod_dataset.RasterXSize
			rows = aod_dataset.RasterYSize
			
			avg_temp_file = "/home/phamha/MetTiff/temp"+str(month)+".tif"
			avg_rh_file ="/home/phamha/MetTiff/rh"+str(month)+".tif"
			avg_preci_file = "/home/phamha/MetTiff/preci"+str(month)+".tif"
			
			avg_temp_dataset =  gdal.Open(avg_temp_file, GA_ReadOnly)
			avg_temp_data = avg_temp_dataset.GetRasterBand(1).ReadAsArray()
			
			avg_rh_dataset =  gdal.Open(avg_rh_file, GA_ReadOnly)
			avg_rh_data = avg_rh_dataset.GetRasterBand(1).ReadAsArray()
			
			avg_preci_dataset =  gdal.Open(avg_preci_file, GA_ReadOnly)
			avg_preci_data = avg_preci_dataset.GetRasterBand(1).ReadAsArray()
				

			pm_datal = numpy.zeros((rows,cols), dtype=numpy.float32)
			for i in range(0,rows):
				for j in range(0,cols):
					if aod_data[i,j] != -9999 and temp_data[i,j] != -32768:
						aod = aod_data[i,j] * 0.00100000004749745
						temp = (temp_data[i,j] + 15000) * 0.00999999977648258
						avg_temp = avg_temp_data[i,j]
						avg_rh = avg_rh_data[i,j]
						avg_preci = avg_preci_data[i,j]
						pm_datal[i,j] = quadratic(aod,temp,avg_temp,avg_rh,avg_preci)
						#pm_datal[i,j] = 25.6965586 * aod + 0.1525928 * temp + (-2.1894375)*avg_temp + (-0.6837351)*avg_rh + 0.4189560 * avg_preci + 70.1740458
					else:
						pm_datal[i,j] = -9999.0
			
			pm_output = os.path.join("/var/www/html/fimo/apom/Product/MYDPM",str(year))
			pm_output = os.path.join(pm_output, aod_filename)
			if not os.path.exists(os.path.dirname(pm_output)):
				os.makedirs(os.path.dirname(pm_output))
			
			gtiff = gdal.GetDriverByName('GTiff')
			output = gtiff.Create(pm_output, cols, rows, 1, gdal.GDT_Float32)
			if output is None:
				print "Output dataset is none"
				exit(1)
			
			srs = osr.SpatialReference()
			srs.ImportFromProj4('+proj=longlat +datum=WGS84')

			output.SetProjection(srs.ExportToWkt())
			output.SetGeoTransform(aod_dataset.GetGeoTransform())
			# output.SetNoDataValue(-9999)
			output.GetRasterBand(1).WriteArray(pm_datal)

			output = None
			
			print pm_output

	except psycopg2.DatabaseError, e:
		print 'Error %s' % e    
		sys.exit(1)
	finally:
		if con:
			con.close()

if __name__ == '__main__':
	main()

